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ABSTRACT 

We use the transfer equation in relativistic form to develop an expansion of the one- 
photon distribution for a medium with constant photon mean free path, e. On carrying 
out appropriate integrations and manipulations, we convert this expansion into one for 
the frequency-integrated intensity. We regroup the terms of the intensity expansion 
according to both the power of e and the angular structure of the various terms and 
then carry out angle integrations to obtain the expansions for the components of the 
stress energy tensor: the radiative energy density, the radiative flux and the pressure 
tensor. In leading order, we recover Thomas' (1930) results for the viscosity tensor and 
his expression for the viscosity coefficient, which are correct for short mean free paths. 
As had been done earlier for the radiative heat equation (Unno and Spiegel 1966), we 
keep at each order in the expansion a dominant portion, but this time one with a richer 
angular structure. Then, after some rearrangement of the various summations in the 
expressions for the moments, we replace the sum of the calculated higher order terms by 
a Pade approximant, or rational approximation (Baker 1975), to provide an improved 
closure approximation for the radiative stress tensor. The resulting radiative viscosity 
tensor may be expressed either as a simple integral operator acting on the Thomas stress 
tensor or as the solution of an inhomogenous, linear partial differential equation. The 
expression obtained for the radiative viscosity tensor applies for media with long, as 
well as short, photon mean free paths. We also develop results applicable for relatively 
smooth flows by using the form of the Thomas stress tensor with generalized trans- 
port coefficients derived by the application of a suitable operator to the bare Thomas 
coefficients. 

Subject headings: radiative transfer; resummation; stress tensor; transport coefficients; 
entropy generation; radiative viscosity; radiative fluid dynamics 



1. Introduction 



The presence of a radiation field will subject an otherwise perfect fluid to dissipation. In the 
limit of short photon mean free paths, dissipation is expressed as diffusion of momentum and energy 
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and is regulated by the transport coefficients of viscosity and conductivity. These were obtained 
by Thomas (1930) who started from the relativistic transfer equation and developed the solutions 
in an expansion in photon mean free path. Since Thomas retained only terms of first order, his 
was a diffusion theory whose results have been extended to include the effects of general relativity 
(Hamecn-Antilla and Suhonen 1971) Thomson scattering (Masaki 1971; Hsieh and Spiegel 1976) 
and Compton scattering (Masaki 1981). 

Unfortunately, a diffusion approximation does not provide an adequate description of the effects 
of radiation on the medium and its dynamics when vigorous radiative fluid dynamics produces large 
fluctuations in physical conditions. In such conditions, transparent regions may form and the photon 
mean free path can become relatively long so that better approximations are needed. One approach 
is to continue using the equations for the first two moments of the radiative intensity, energy density 
and flux, and closing the system by providing an improved approximation for the pressure tensor. 
A step in this direction is to use the Eddington factor methods (Mihalas and Mihalas 1984) as has 
been done by a sequence of authors (Castor 1972; Fu 1987; Minerbo 1983; Dominguez 1995, 1996, 
1997). As Castor (1972; see also (Mihalas 1983)) has reported, the scheme seems to be usable 
for long photon mean free paths, but it does call for the introduction of additional information. 
Typically the variable Eddington factor is to be found from numerical solutions of the transfer 
equation in related conditions. It is therefore desirable to seek a self-contained approximation for 
the pressure tensor and we shall provide one here. 

We shall here use an approach that originates, in part, with a method used for closure in the 
study of thermal effects of radiation (Unno and Spiegel 1966). In the treatment of the thermal effects 
alone, the Eddington approximation gives rather good results. This is so because the Eddington 
closure relation for the radiative stress-energy tensor is the result of resummation of terms of all 
orders in an expansion in photon mean free path (Unno and Spiegel 1966), and it captures aspects 
of both optically thin and thick structures. When we use the Eddington approximation we get a 
radiative heat equation (for a static medium) that is schematically of the form (Unno and Spiegel 
1966): 

d t s = Kl vv 2 s, v= (i- K2 v 2 ) -1 , (i) 

where k\ and K2 are suitable "diffusion" coefficients. The Laplacian, which is the usual diffusion 
operator, is here augmented by the additional smoothing operator V. Of course, this is the bare 
essential of such a result, whose simplicity depends on a number of idealizations, most notably that 
the mean free path of photons is a constant, and on the neglect of retardation effects. Nevertheless, 
T> carries much of the qualitative character of the radiative smoothing process. If we develop the 
smoothing operator in equation (1) for large scales (small V 2 , qualitatively speaking) we obtain 
W 2 ss V 2 + K2V 4 . This illustrates why we may recover the full expression for V from a Pade 
approximation based on only the first two terms in the expansion. 

On the other hand, for the case of viscous stresses, the Eddington closure approximation fails 
to produce a physically acceptable result (Anderson and Spiegel 1972) and an improved version 
is needed. The aim of the present paper is to devise a closure for the second order moments of 
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the radiation field that may be used to deal with both the dynamical and thermal effects of the 
radiation field for both optically thick and thin structures. Similar techniques have recently been 
used in kinetic theory where an improved viscous stress tensor has been obtained by resumming 
terms up to only second order in an expansion in mean free path (Rosenau 1989, 1993; Slemrod 
1997). The resummed forms, when expanded, contain terms of all orders, as in the earlier radiative 
thermal study. 

Starting from the manifestly covariant form of the transfer equation for a moving medium, we 
develop the solution of the transfer equation in photon mean free path. Then we group terms within 
this expansion according to angular complexity. This leads to multiple sums and we rearrange the 
series so as to be able to compare terms of like character. From each order, we take the dominant 
term in power of the mean free path. Then we choose the leading portion at each order in terms 
of angular complexity, in what is really a Galerkin truncation. The resulting infinite sum is then 
converted to finite form by a resummation, or Pade, technique (Baker 1975) to provide a compact 
form for the radiative viscosity tensor. 



2. Basic Formulas 

2.1. Transfer Equation 

To include the effects of differential motion of a medium on the transfer of radiation, Thomas 
(1930) used the transfer equation in relativistic form. In this way, he readily included effects such as 
aberration and Doppler shift and was able to obtain the correct radiative viscosity where others had 
failed (Mihalas 1983). Even if one is not concerned with relativistic fluid velocities, it is convenient 
adopt the relativistic formulation and also to follow Hazelhurst and Sargent (1959) in using the 
covariant form of the equation of radiative transfer for the case of moving media. This approach is 
especially effective if one makes use of the one-particle distribution function for photons, f(p,x), 
where p (or p^) is the four-momentum of a photon and x (or x^) is its location in space-time. 
According to the way one proceeds, / is a scalar (as here) or a scalar density. 

If we neglect gravitational bending of light and refraction by the medium, we may write the 
transfer equation in the general form (Anderson and Spiegel 1972) 

p^f = p(a-f3f) (2) 

where stands for d/dx^ 1 and summation over any repeated index is understood. Here a is the 
rate of injection of photons with momentum p by emission and j3f is the rate at which photons 
are removed by absorption; both rates are per unit mass. The quantity p is the number density of 
absorbers times the rest mass per absorber. 

Unlike /, the specific intensity, 2, is not a scalar. Its transformation properties, as given by 



-4 - 



Thomas (1930), may be deduced from the definition 

^ h^ 3 f 



c 3 



(3) 



where h is Planck's constant, c is the speed of light and v is the frequency. We shall use units in 
which c and h are unity. 

The transformation properties of the other quantities in the usual transfer equation are ob- 
tained by comparing that equation in the local rest frame of the matter with equation (2). For 
our purposes, there are two important reference frames, a basic inertial frame (such as that of a 
stationary star) and a frame locally comoving with the matter. The four-velocity of the matter, 
u^, has components (1,0,0,0) in the locally comoving frame. 

The invariant 

v = u% . (4) 

is the photon frequency in the local rest frame of the matter. By comparing the transfer equation 
in a general frame with that in the comoving frame, Thomas (1930) found (3 = vk(u), where k is 
the opacity, and a = i> _2 j(z>), where j is the emissivity. Since a, (3 and v are scalars, we then know 
the transformation rules for j and k once we introduce the (relativistic) Doppler relation between 
general frequency v and v. 



2.2. Radiative Stress Tensor 

A basic object of radiative fluid dynamics is the radiative stress tensor, T^ v . It may be defined 
as an integral of the distribution function over momentum space: 

TV = J p^fdP , (5) 

where dP is the covariant volume element in momentum space. With v as the radial coordinate 
in spherical polar coordinates in momentum space the element of volume in momentum space is 
dP = vdvdVL, where dQ. is the element of solid angle (Landau and Lifshitz 1984). We may then 
factor the frequency from and define the null vector 

n^=p^/D, (6) 

so that the stress tensor is 

= / d£ln> x n v l (7) 
where we have introduced the frequency-integrated specific intensity 

1= / v^fdv . (8) 
Jo 



-5 - 



The projection operator 

h v» = v ^ _ u » u » ( 9 ) 

picks out space-like components in the frame of the matter, where rf v is the Minkowski metric, 
diag(l, —1, —1, —1). Then we may introduce the spacelike unit four-vector, 

l» = h ^n u , (10) 

so that 

= + V . (11) 

We see that V 1 must be perpendicular to and that = —1. Then with the introduction of 
(11) we arrive at a standard way to write the stress tensor, namely 

T» v = P» v + F^u" + u^F v + Eu^u u (12) 

with 

e= idn , f^= ridn , = / m v idn , (13) 

where these integrals are over all solid angle. 

By taking moments of the transfer equation itself, we derive equations for moments like those 
in (12). This procedure leads to an infinite hierarchy of moment equations that may be truncated 
by introducing a relation amongst the moments, E, F^ 1 and P^ v . Such closure relations can be 
found by assuming some particular forms for the angular dependence of I, usually expressed in 
terms of the V 1 or n M and this requires angular integrals of their products. Though we do not 
follow this procedure, we do need such integrals and so we give a list of formulas for the angular 
integrals in Appendix A. In the simplest case, when the isotropic approximation is introduced for 
the intensity, we need the integral 

n u dn = -^-h^ (14) 

for the evaluation of P^ u that is used in the Eddington closure approximation (Anderson and 
Spiegel 1972). 



2.3. The Matter and the Total Stress Tensors 

To describe the matter macroscopically, we adopt, with Weinberg (1971), the procedures given 
by Eckart (1940), who noted that a vector may be written as = {u a V a )u^ + haV a , as in 
equation (11). A similar decomposition is possible for tensors and, in particular, for the radiative 
stress tensor, we see that 



E = u a u p T a P , F^ = u a T a ^h^ , P^ = h£h"p T aP . 



(15) 
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As in Eckart (1940) 's approach, we may use the form of (15) to define the three important 
moments of kinetic theory in terms of the matter stress tensor. If we denote the stress-energy 
tensor of the matter as T^, the mass-energy density, the matter flux and the pressure tensor are 
given (as in (15)) by 

£ = u a u p T a ^ J* = u a r al3 h^, = Kh%T a P . (16) 

To concentrate on radiative effects, we assume that, in the absence of interaction with radiation, 
the matter would be a perfect fluid, for which = and V^ v = Vh^, where V is the partial 
pressure of the matter. The matter stress tensor in that case would be 

T"" = Eu»u v - VW V . (17) 
The total stress tensor (matter plus radiation) is 

and the equations of motion of the mixture are expressed as 

9^ = 0. (19) 



3. Diffusion Theory 

For a discussion of the extension of Thomas' theory, a brief restatement of his results is useful 
for bringing out some of the issues to be faced. To provide one here, we introduce the source 
function, S = a/[3, and the local photon mean free path, e = (pk)~ 1 where k is generally a 
function of v as well as of the state variables of the medium. The mean free path appears explicitly 
in the transfer equation (2), which we rewrite as 

i = S-en»f# (20) 

where is defined in (6). This way of writing the transfer equation is meant to take advantage 
of the smallness of e, which is assumed in Thomas' theory. The fact that e appears in front of the 
derivative means that the perturbation theory based on its smallness will be singular. Generally 
speaking, singular perturbations expansions give rise to (divergent) asymptotic series but these may 
be made useful by appropriate (re)summation techniques (Bender 1978) such as we shall use here. 

In first approximation, we have / w S. If we iterate on this once in (20) we get Thomas' 
approximation, 

/ = S - en^ . (21) 

To compute the stress tensor, we first perform the integration over frequency after multiplying 
by z> 3 to obtain an approximation for the frequency-integrated intensity, 

poo 

/ = < S-n^ / iPeS^dv, (22) 
Jo 
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where we have introduced the integrated source function 

S(T) = [ D 3 S(i>,T)dD . (23) 
Jo 

The integral in formula (22) is a simple example of the kind of integrals that the higher order 
development entails, so it is worth discussing its evaluation. 

In Thomas theory, D 3 S is the Planck distribution and it depends explicitly on v and temper- 
ature, T, from which it derives an implicit dependence on position and time. Therefore we may 
write 

S tli = v^doS + T^d T S . (24) 

On recalling that v = and that does not depend on coordinates, we obtain v jfl = p v u v ^. 

Then, since u v u v ^ = 0, we have that 

z> M = vn v u v il . (25) 

If, as in local thermodynamic equilibrium, the source function depends on v and T only through 
the combination 9/T, we have 

vdo = -Td T (26) 

and so (24) can be rewritten as 

S,n = (Tfj, - Tn'uw) 8 T S . (27) 

The integral in (22) then becomes 

POO f'OO 

/ v*eS^ dD = (Tp - Tn v u v ^) / D 3 ed T S dD . (28) 
Jo Jo 

If we pull the T derivative out of the integral, we get an extra term involving the T derivative of 
e. On the other hand, we can simply introduce a weighted frequency integral through the relation 



\ iPs, T {v,T)edv . (29) 
Jo 



The quantity e is a mean of the inverse of the absorption coefficient with respect to 5 t in the spirit 
of the Rosseland mean absorption coefficient. Then we find that, for TS^t = 4S 1 , 

I = S- enPS^ + AenVSu^ . (30) 



In the higher order theory to be discussed below, we can use similar tricks by introducing a 
variety of mean absorption coefficients, as is done sometimes in the theory of stellar atmospheres. 
More often than not, all these are assumed to be equal in the final analysis. To avoid such purely 
formal maneuvers, we shall assume in the presentation of higher order theory that the medium is 
gray. 
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Then, to this approximation, partial integrations such as we have just described give us the 
known forms of moments of the intensity in the Thomas approximation: 



E = AttS - Aire 
4ir 

F M = — eh» a 
3 



u a S, a + ±Su% 



(31) 

S, a -ASu^u a ,p , (32) 



1 1 £i 

pi* = --hTE + -^-eSu p , a T^ pa , (33) 
o ID 

where we have used the relevant formulae from Appendix A and have introduced 

3 

From the expression for the flux we see that, since S is proportional to T 4 , there is the usual 
conduction term proportional to the (four) gradient of T with a conductivity coefficient 16neS/ (3T). 
This is the diffusion limit and it is much less good than the Eddington approximation, for which 
the flux is proportional to the gradient of E. With the Eddington approximation, thermal times 
are well approximated over all ranges of optical thickness (Unno and Spiegel 1966). The second 
term in (32) has been considered to be a purely relativistic term (Mihalas 1983), but a classical 
analogue exists in the kinetic theory of rarefied gases (Chen PhD thesis; Chen et al. 2000). 

In the expression for the pressure tensor, the first term, to the order given, does mimic the 
Eddington approximation and we have in addition a viscous term. From the standard form of 
the pressure tensor, we can read off the viscosity coefficient given by Thomas, 16e7rS'/15. Like 
the Thomas expression for the flux, this is also a diffusive approximation good only for short 
mean free paths. Our aim here is to improve on the non-diagonal viscous tensor in the same way 
that the Eddington approximation improves on the conduction term. In a subsequent paper, we 
shall describe an improved approximation for the diagonal terms for the case where scattering is 
important. 



4. The Expansion of / 

To go beyond the Thomas approximation, we could continue the iteration process described 
in the previous section, but we find it more effective to proceed by formal expansions. Our small 
parameter will be a typical value of e, call it e. Then we form the expansion 



oo 
m=0 



(35) 
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which we introduce into the transfer equation. On demanding that the expanded equation is 
satisfied term by term and prescribing that e is everywhere of order e, we get 

/(0)=«5, (36) 
as in comoving local thermodynamic equilibrium and, for the higher orders, we obtain 

/(m+1) = --^/(mV • (37) 

In the first order, as in diffusion theory, the complication caused by having e in the formulas is 
slight, but in higher orders, it makes the development very cumbersome especially on account of the 
presence of derivatives of e in the expressions of fr m ) • These may be converted to derivatives with 
respect to rest frequency and state variables as for the source function but still, in the performance 
of the integrals over frequency, it soon becomes necessary to introduce a plethora of mean absorption 
coefficients, and even tensorial mean absorption coefficients (Anderson and Spiegel 1972). Though 
we do not report on this here, we have retained variation of e in calcualtions out to the third order 
in e in (Chen PhD thesis). This is far enough for one aspect of our approach to be carried out, 
though we shall not go into this here. 

More generally, it is possible to make allowances for the dependence of e on by a formal 
device like used in stellar atmosphere theory. We may introduce new coordinates through 
the transformation edy^ = edx^ where e is a mean (over frequency) absorption coefficient. This 
formal simplification is analogous to a transformation to the local rest frame of the matter. But 
since both e and the u M are dependent on the radiation field, such implicit definitions of coordinates 
necessitate, at the end of the calculations, some rather hard inversions. Moreover, we would need to 
take frequency means to use such transformations effectively. We shall avoid all these complications 
here and will stick to inertial coordinates, requiring that e is constant. In that case, we obtain 

/(m+l) = - n>1 f(m),n i (38) 

or, in terms of <S, 

f. > — (--\\ m 

Ds m ' 



f(m) = (-ir-^, (39) 



where 

— = n»d» . (40) 

To perform the indicated differentiations, it is useful to have a formula for the derivative of 
the null vector, n M . On differentiating = n^v and making use of (25), we find that 

n* v = -n^-nPup^ . (41) 

We can then develop explicit formulas for the /( m ) in terms of the source function and these together 
lead to a series providing / as a functional of S. This is tantamount to what Chandrasekhar (1960) 
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has called the formal solution of the transfer equation. (Its integral expression in terms of optical 
distance is helpful in treating the case of variable e.) The idea of the formal solution is that the 
source function and the absorption coefficient are material properties of the medium and are to 
be considered as given for the purposes of solving the transfer problem, though the situation is 
reversed when we try to solve equations for the properties of the medium. So we need to decide 
how the source function is to be specified. 

In the simplest case mathematically, we would be given S explicitly as a function of space and 
time. From this, we would compute frequency integrals to get the integrated intensity through 
formula (8). But it is more representative of the general situation to assume, as we did for diffusion 
theory, that <S is expressible in terms of T and v. Though the source function does not depend on 
direction, we do have effects of angular distribution entering in through the transfer equation itself. 
To make allowances for that also, we write 

= Tfj, 8 T + P p u p ^ do + rfjn* (42) 

and we are led to a generalization of the differentiation formula implicit in (30): 

8^ = (T^ - TnP Up ^) 8 T - rftfu^ d n „ . (43) 

On noting that the quantity /( ) = S does not depend on direction (n^) we can evaluate f^y 
We find using (43) that 

= -StTX + TS tT u v ^n u . (44) 

If we use this formula to compute the moments, we recover the results of the previous section with 
e in place of e. 

From (39), we see on using (41), that 

/ (2) = n»n v S^ v - n»n v n<>u p ^S, u . (45) 

We then have to cope with the higher derivatives of the source function S. If we continue to assume 
that S depends only on i>/T, since S does not depend on n^ 1 , (43) gives us 

S, li = TS, T [(lnT) tll -{lni>\ li ]. (46) 

A similar, but lengthier, calculation leads to 

S tlu , = TS, T {[{lnT)^-(lni>)J(lnT)^ + (lnT)^-(lni>)^} (47) 
+ [T 2 SM^T), V - (TS, T + T 2 S, TT )(lnD), u ] [(InT)^ - (lni>)J 

where we have used (25) and (41). 

We then find that (45) takes the form 

/ (2) = - n»n v rtjX$ + n»n"nOn°TW c (48) 
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with 

F p v ] = s ,tT,hu + S t TrT^T v 

FgSo = (3T5 r + r 2 5 TT )n^n^. (49) 

This way of writing // 2 ) is useful in computing terms in the next order and from (39) and (40) 
and the forms of /(i) and /( 2 ), we may surmise that 

2m 

f(m) = J2(- 1 ) k 4T-% k n^---n^ (50) 

k=m 

where the quantities T are functions of S and its derivatives with respect to T; that is, they are 
functionals of S. (If we had not assumed constant photon mean free path, they would be functionals 
of loge as well.) The calculation of the distribution function at any order then boils down to the 
calculation of these functionals, which is a straightforward, if demanding, task. For m = 3, we are 
led to the results 

Ffvl = ^,TTTT^T v T a + 3S iT TT pu T a + S :T T ifiaL , (51) 

•^S = TS, T u a ^ v + (4«S T + 3TS yTT )u a ^T u + (6<S T + 3TS, TT )u a ^T pu 

+ (85 T t + STS^TT^pT^Tv (52) 

^Sp = u a A(WTS, T + 3T 2 S, TT )u^ p 

+{15S, T + 18TS, TT + 3T 2 S,ttt)u^T p ] (53) 

= [MTS t T + 9T 2 S, TT + T 3 S, TTT ] U ,, u u p>aUa>p (54) 
We thus have all the terms up to third order. 

In this way we see that we may order the terms by increasing angular complexity. In that case, 
if we go up to four factors of n^, we have an expansion of / in the form 

-e 2 + ejffi*, + e 2 ^t%]n^n^ + ■■■. (55) 

We may compute the coefficients T at any order using this procedure. 



5. The Intensity Expansion 

An intermediate step on the way to computing the moments that make up the radiative stress 
tensor, is the computation of the (frequency) integrated intensity defined in (8). Its expansion is 

oo 

/ = £ e m I (m) , (56) 

m=0 
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and when we introduce (50) into (8) we obtain for the coefficients in (56), 



2m 

^ ) = E(- 1 )%™^ 1 ---^ 

k=m 



where 



Y, 



(m,k) 



(57) 
(58) 



With the expressions for the J 7 ^)^ given in the previous section, we may develop explicit 
expressions for Y^.' k J k . We do this by working out a recursion formula for these coefficients. First, 
we introduce the expression (50) for fr m ) into (38) to obtain 



2m 



/ (m+ D = -^E(- 1 ) A 



k=m 



-r( mfc ) r,Mi . . . i -p(mk) , M1 _ \x h \ 



Since from (41), we have 

(n w ■••n' ifc ) )At = -kn p u p ^{n^ ■■■n fMk ) , 
the only real complication comes from the derivatives of the J- . 



(59) 



(60) 



Next, we reexamine the derivation of the explicit formulae for the and see that their 
dependence on position and time comes about through their dependence on z>, on T and on u^, as 
well as the derivatives of these quantities with respect to So each T depends on v and on N 
other variables that we shall call V A , where A = l,2,...,N. The V A are T ;(U , T^, u^, and so 
on. Let 

A<™>*> = V* dy A 

where summation over repeated A is understood. Then 



On recalling that = Dn p u p ^ we can rewrite (59) as 

2m 



(61) 
(62) 



f {m+1) = E(- 1 )' C+lpnPn %.M^ 



(mk) 



/'I 



k=m 
2m 



+ ^{-l) k+1 n»K^ -knPn^Up^ 



(63) 



k=m 



To compute the /( m ), we carry out the frequency integral indicated in (8) and integrate by 
parts in the term with <% to obtain 



2m 



k=m 



W) = E(- 1 )' C (H4)nV Mw -/A^) y^n- 



(64) 
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Now v is no longer in the problem, and so (62) simplifies to A^ 1 '^ = d^; we shall make this 
replacement in the rest of the development. 

Formula (64) has two terms, each having a different angular structure. To keep to our general 
attempt to order things by angular structure as well as by mean free path, we make a further 
rearrangement. So far, in the expressions for the for different pairs (ink), we have had need only 
of values for which m < k < 2m. Therefore we are at liberty to impose that 

•^Mr-jL = for k < m and k > 2m . (65) 

Accordingly 

Y^: k J k =0 for k < m and k > 2m . (66) 

Moreover, for either m < or k < 0, we define Y^™.' k J k to be zero. 

Now we can adjust the indices in the summation (64) by replacing k in the first term with 
k — 2 and k in the second term with k — 1. We then have 



2(m+l) 
\ ^ I i\k 



W) = E (-!)* [( k + 2)4 m 42 ^ M + YfctS ,» k \ n^-.. n "> (67) 

k=m+l 

and, on comparing with (57), we observe that 

v (m+l,k) _ v (m,k-l) ^ v (m,k-2) , fi o\ 

which will be used in the sequel. The Y sequence is to be built up from As we see from 

(58), this is the integrated source function, that is 

y(°>°) = S. (69) 

Hence, as seen above, / is a functional of S as well as of and the angular variables. 

6. The Viscosity Tensor 

Once we have calculated the integrated intensity, as in the previous section, we can evaluate 
the physically important angular moments, 

CO OO CO 

m=0 m=0 m=0 

Since the Y defined in (58) are independent of angle, the terms in these expansions are 

2m „ 

E (m) = E (- 1 ) fey Si J ™ V • • • «"* dn ( 71 ) 

k=m 
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F (m) = E / l " nV2 ---^ ( 72 ) 

fc=m 

2m „ 

P ( m) = E (- 1 ) fey iSl y i^"" 1 "" 2 ■ • • dn . (73) 

k=m 

Though we shall quote some results about the approximations to E and F p , our main interest 
is in providing an improved treatment of the viscous effects for transparent regions of the system 
under study. So we shall detail only the calculation of an improved approximation for the pressure 
tensor. If we introduce the Eddington closure for the diagonal terms in the stress tensor, we may 
write 

pv» = — Eh^ + E^ , (74) 
3 

where E pv is the viscous stress tensor, whose expression the main goal of this work. As mentioned, 
in a later paper, we shall show how the first term in (74) can be improved considerably. 

To seek a representation for E pu that is suitable for both long and short photon mean free 
paths, we begin by expanding the viscosity tensor as 

oo 

= E • ( 75 ) 

m=0 

From (71) and (73) we find the coefficients in this expansion and so can write 

2m 

S (m) = E (-l) k Y^: k J k f2^-^, (76) 

k=m 

where 

QUVM-Hk = / ^ x v + _^) n Mi n M2 . . . n W= dn ( 77 ) 

Since the trace of IH V + is zero, and E pu , are traceless, symmetric tensors. 

For the Ys we have the recursion formula (68) given at the end of the previous section, while 
to reexpress the angular integrals usefully, we can refer to the formula of Appendix A. We see from 
these that 

W = 0, n^P = and W pa = — T ^ pa , (78) 

15 

where T pvpu is defined in (34). To extend the results of Thomas we next generalize r^ p<T into a 
series of orthogonal polynomials that represent the expansion of f2. 

If we use (A6), we can develop the products of the n M and express fi/"W"Mfc i n terms of the 
moments of P. Then we have 
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where 

^Yftvm...fi t _ ^iwn 1 ...n l _| y^v j^fn...^ ^ggj 

3 



and 



J l^F 2 ■ ■ ■ l^dQ . (81) 



The integral for an odd number of factors is zero so, for odd k, M fM1 '"^ k = 0, and hence 
V^/ii^i-.-Ws = 0. Therefore we may rewrite (79) in a way that will be convenient later, namely, 

f-1 

2 /2?'\ 
i=o ^ ^ ' 

where [|] means the integer part of f . 

When k = 1, the upper limit of the sum is % = 0, and = W^u 9 ' . Because M = 4ir and, 
according to (14), Af" = -%hT, we have Ty Ml/ = and ri MJyp = as just noted. In the same way, 
at the next level, from (A2) we see that 

which is equivalent to the third expression in (78). We then define a generalized r as 

3 

To introduce this into (82) we rewrite (84) as 

where it is understood that W viiq is zero if either p or g is negative. Hence, for £ = 1 we recover 
(83) and we may then proceed recursively to express W in terms of r. When £ = 2 we have 

3 

From this we see that 

^■jUI//il-/X2p _ ^ I __ I T tU'Hl—tl2qf l H2q+ltl2q+2 . . . fo^p-l^p _ (g^ 

When we put (87) into (82) we obtain the desired expansion of O in terms of the generalized r: 

= E( 2 ;)£(-r" 

i=0 v 7 j=0 v 

x T IU>(lJ,l—H2j . . . ^2i-lM2i u ^2i+l . . . ^Mfc) _ (gg) 

The generalized r are useful for developing the angular structure of tensors and we find them more 
suited to this problem than the spherical harmonics of the usual moment schemes. 



7. Ordering the Summations 



In a calculation that is in some ways the predecessor of this one (Unno and Spiegel 1966), 
a series of angular moments was rearranged and resummed to provide a closure approximation 
for the thermal terms. These terms lacked off-diagonal elements and so gave only the Eddington 
approximation, which is nevertheless a great improvement on the diffusion approximation. The 
moral is that it is very effective to ensure that all orders in mean free path made themselves felt in 
the result. However, since the approximation used in that earlier work fails to capture the viscous 
effects, we here propose one that retains the off-diagonal terms and leads to results good both for 
long and short photon mean free paths. To bring out the nature of the rearrangement we propose, 
let us shift the index k in (76) to k + m so that (75)-(76) becomes 

oo m 

SF = EE T t) ( 89 ) 

m=0 k=0 

where 

T K*) = {-l) k+m e m Y^ m ^l^--^ . (90) 

When we write out the T^ fe as an array of tensors, with m designating the rows and k 
designating the columns, it takes the form 



1 (0,0) 








1 (1,0) 







1 (2,0) 


X (2,l) 


1 (2,2) 



(91) 



Thus, the Tf u , , form a lower triangular matrix of tensors. 

' (m,fe) ° 

The first sum in (89) is taken over the m th row in the array of Ts and the second sum is over 
the columns. In the spirit of the earlier derivation of the Eddington approximation by resummation, 
we rearrange the terms so that the sum is first taken over columns, then over the rows. That is, 
we can sum first over the m. Since, for fixed k, the sum over m begins with the value m = k, it is 
then convenient to shift the index m by k. This simple rearrangement amounts to 

oo m oo oo 

V y t' o/ r?" . (92) 

/ j / j (m,k) / j / j (m+k,k) v ' 

m=0 k=0 k=0 m=0 

This device, which we shall employ again, turns (89) into 

oo oo 

^ = EE T Um ( 93 ) 

fc=0 m=0 
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so that the expansion of the viscosity tensor is now 

oo 

^-E efes E ( 94 ) 



fc=o 

with 



oo 



_ V (-If e ro y' m+,i ' m+2 ''(]Wi ,, Ma (95) 



m=0 

The two expansions for H^, (75) and (94) look similar but, when we compare them carefully, 
we detect an important difference: only the series for Hj^j contains terms of all orders in e. We 
see also that (93) involves fc w which in turn requires that we be more explicit about the 

Y^^^ltk^ m (95)- The relevant development can proceed with the help of the iteration formula 
(68). We first use (66) in (68) with m — > m — 1 and k = m to obtain 

v (m,m) _ v (m-l,m-l) , Qf ,-. 

Iteration of this formula gives 

v (rn,rn) _ v (0,0) fQ7 , 

More generally, for k > 0, we may use (68) to find 

^mi-'-mL = (2^ + 2)^[ 1 .../i2 fc _ 2 ^ u ii 2 k-i,^2k ■ (98) 

Again, using (68), we see that 

v (fc+l,2fc+l) _ v (fc,2fc) o^(fe,2fe-i) fqq ^ 

and so on. For the special case of k = 0, we obtain (97). 

Now we can use the rearrangement trick again by noting that if we put (98)- (99) into (95), it 
takes the form 

oo rrt 

S W = Yl J2 X (m+k,k,i) ( 10 °) 
m=0 i=0 

where the expressions for the X can be obtained by comparing to (95). We do not need their 
explicit form since we need only to observe that by a repetition of the formula (92) we get 

oo oo 

-[k] =EE X (m+k+i,k,i) ■ ( 101 ) 

We may then carry out the necessary substitutions, and separate out the exceptional case k = 0, 
to obtain 

oo oo oo 

= ^2 E E ( — l) m+i e m+k+i (2k + i + 2)S7 mil '" M2fe + m +^ 

k=l i=0 m=0 

. (fc+i-l,2fc+i-2) 

1 1 Hl---H2k+i-2 a pt2k+i~\,^2k+ 



o 

+ ^(-l^^y^QW"^ (102) 



,M2fc+l + i-"M2fc + m + i 

.(0,0) 

m=0 
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This completes the formal expansion procedure for the viscosity tensor, which we can now 
express in terms of the generalized r and the coefficients of the intensity expansion. For both of 
these sets of quantities we have recursion formulae. We next begin to make approximations to 
prepare the way for representing the results in terms of rational approximations to the series. 



8. Approximating the Orders 

So far we have been proceeding with formal expansions in mean free path and making conve- 
nient rearrangements of the series obtained. The result is that terms in the series depend on e so 
that terms of similar character can be compared side by side. Now we turn to some approximations 
to make the series more manageable. What we shall propose in this section is the neglect of those 
terms in each order of the formal expansion that are clearly (and asymptotically) smaller than 
their fellows in the same order, as in the earlier study of the radiative heat equation. There it 
was observed that by retaining a dominant term in every order one is led to a useful outcome (in 
that case, the Eddington approximation). That approach is what we mean by 'approximating the 
orders.' 

We note that in (102) we may change the order of summation, since the sums are each taken 
to infinity in that expression. So we may rewrite (102) more compactly as 

OO s OO OO N 

= E(- £ ) m EB- 1 )'^) + Y$l% m W^-^ (103) 

m=0 ^fc=li=0 ' 



where 

A Zi ki) = ( 2k + i + 2)n^ 1 "^ 2k +^+ i 



Y (k+i-l,2k+i-2) 

1 Hl---H2k+i-2 U M2fc+i-l>M2fc+i 



(104) 



We operate first on the double sum of the summand of (103) by separating out its leading 
term for k = 1 and i = 0. The double sum can then be rewritten as 

OO OO OO OO 

^,1,0) + E Ei- 1 ) 1 ^ 1 ^) + E efc ^,o) + E(-^ i+lA ^,i) ■ ( 105 ) 

k=2 i=l k=2 i=l 

In the double sum of (105), we shift the indices so that k becomes k + 2 and i goes over into i + 
Similar shifts can be carried out for the last two terms of sum. Then (105) becomes 

OO OO OO 

^,i,o) - ^ E Et- 1 )^:,^) + e2 E e " fe +2 ,o) - (-1)^,^1)] • doe) 

fc=0 i=0 k=0 

Since the last two terms in this expression are down by at least a factor of e from the first term, 
from which they are otherwise not dissimilar, we may neglect them. 
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Next, let us look at the second term in the summand of (103). A simplification of this term 
can be based on the fact, seen on inspection, that W / = 0, together with the result given in (78) 
that = 0. With the first two terms gone from this second sum over m it becomes 



J2 (-e) m Y^ m n^-^ = J2 {-€) m Y^^Jl^-^ . (107) 



'(0,0) Q/iiz/n— fi m _ ( _^\m v (0,0) 

m=0 m=2 

By shifting m to m + 2, we can turn (107) into 

oo 

^ ( _ er+ 2y(0,0)^ +2 ^ 1 ... Atm+2 ; (10g) 
m=0 

and so, we can write (103) as 



OO s OO OO 

m=0 k=0 i=0 



m,k+2,i+l) 



+ e 2 Y^% m+2 n^-^ + ^ . (109) 

As we have stated, we are going to neglect the middle term in this summand and now we see that 
we may also neglect the last term. Our expression for the radiative viscous stress tensor may then 
be approximated as 

oo 

•ET = ±eJ2 {-e) m WW"»» [^ (0 ' 0) <wU...^ (HO) 

where each of the neglected terms is O(e) compared to its counterparts. We note that the term for 
m = of this expansion is the Thomas approximation to the stress tensor. 



9. The Viscous Smoothing Operator 

The thermal smoothing operator in (1) is a rational fraction of polynomials in the Laplacian 
operator. It was obtained by first expanding the original integral operator of radiative smoothing 
in a series of powers of photon mean free path and approximating each term according to various 
criteria, especially simplicity of geometrical structure. In attempting a similar approach here for 
the viscosity problem we encounter a more complicated development. We have grouped terms in 
our expansion by a combination of order in the mean free path expansion and geometric complexity. 
This approach has permitted us to compare terms of like character and then to keep in each cluster 
of terms the leading ones in powers of mean free path. Now we wish to examine the remaining 
terms and keep only the geometrically simplest terms, as was done in the derivation of T> of (1). 
We shall also make some further approximations in terms of e. 

As always, much depends on y(°>°) = S (see (69)), which we now introduce into (110). We 
then observe that the viscous stress tensor takes the simple form 



(111) 
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where the smoothing operator, ID, the centerpiece of this investigation, is given as 

oo 

jQHupc = ^(.^flW*^ ^ . ( 112 ) 



m=0 



We now wish to extract the geometrically simplest part from each term in this expansion. To do 
this, we recall formula (88) for fiW'^j in which the leading nonvanishing term in the series 
involves the tensor T^ vpu that also appears in the standard form of the viscous stress tensor. The 
origin of the general r has been discussed in the derivation of (87). 

To obtain the usual form of the Eddington approximation, one normally omits higher order 
angular moments of the intensity field in series such as we have obtained. But the conventionally 
neglected moments contain projections onto r^ upa and so that form of truncation misrepresents 
the viscosity tensor. That is why we developed the r tensors — by retaining only the leading 
angular approximation in these quantities we keep just the angular quantity that is central to 
the representation of viscosity. Therefore, we here select from (88) only the geometrically most 
elementary part of the term in each order and write 

[~] £ 1 

^nvfii-fim _ ( _ ~) T^ u ^ lP2 h ps,pA ■ ■ ■ h^ 21 - 1 ^ ■ u^- 1 ■ ■ ■ u Pm ^ (H3) 

In allowing contributions to each term from only the leading generalized r we have made what may 
be called a Galerkin truncation. 

Next we make an index shift j — > i — 1, so that (113) becomes 

j- m — 2 j 

QUVfJU-Vm — /2j + 2\ / r ^l/Oi^2/^3^4 . . . ^M2j + 1 A*2j+2 u /*2j+3 . . . ^m) (114) 

fr, v "' / v v 

This in turn can be written as 

It! j 

QHl/pam-lim _ ST- I 2 j + 2 \ ixvl^pa-, M1/i2 _ _ _ ,H2j-lH2j u H2j+l . . . /m) ("115") 

~ \m + 2j V 3j 1 > 

With this approximation for f2, the viscous smoothing operator is 

rial 

[ 2 J /o„_i_o\ / ^\0 



= E(-rE( 2 i: 2 2 )H) 



m=0 j=0 

x T i*>QxT h nii» . . . ^-iwj^wj+i . . . U ^)^ i ... /Jm . (H6) 

Because 9 Atl ... Mm is symmetric, the symmetrization of the /ij is redundant and may be omitted. 
However, this is not the case for indices p and a. If we approximate (116) by keeping only the term 
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with factor T^ upCT after symmetrization, there is an extra numerical factor ( m+2 ) 2 (m+i) • Taking all 
these into account and interchanging indices j and m in the summation, we obtain 

j=0 m=2j v ' x 

where 

lL {j , m} = h^ 2 ■ ■ ■ h^-^u^+i ■ ■ ■ ^+ 2 ^ M1 ... Mm+2j . (118) 
When we replace the index m by m + 2j (117) becomes 



oo oo 



3 



JD^/XT = T ^P<r Y Y £ r ( 2j + 2 J (-e) m+2 J (--)' iLr,- m> . (119) 

We next want to simplify the differential operator 1L by approximating it with 

Direct evaluation shows that -2^{o,o} = ^{0,0} ; -^{1,0} = ^{1,0} an d -^{0,1} = ^{0,1} j an d this means 
that we would be making approximations only in the higher order terms in (119) in replacing 1L 
by C That is, this replacement leads to the neglect of terms of relative order e 2 in only the higher 
terms, as we detail in Appendix B. 

The advantage in making the approximate replacement of JL is that we may write 

% m} = (h^d af3 y(u a d a r ; (121) 

this lends itself to the rational approximations that we wish to develop for the viscous smoothing 
operator, which we may now write as 



00 / ,2 \3 00 9 / 2i + 2 \ 

jD^pa = I --h a ^d a0 ) V- J )(-eu a d a ) m . (122) 



10. Rational Approximation for JD^^ 

An asymptotic expansion such as we have developed can be made more useful by the technique 
of resummation (Baker 1975; Bender 1978), of which there are various forms. The general idea is 
to replace the series by a rational function of operators in one of several well documented ways. 
For a given function g(e) with series representation 



00 

i=0 



(123) 
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the Pade approximation of order [m, n] is written as the rational function 

*lm,n]W)] = ^ (124) 

where P m (e) and Q n (e) are polynomials in e of degrees m and n, respectively. The coefficients in 
these polynomials are to be fixed according to the condition that the first m + n coefficients of the 
expansion of 3f?[ m n ] should match the first m + n coefficients in the series for g{e). 

This procedure, for all its wide acceptance, has not been rigorously justified, but a considerable 
experience shows that most of the reasonable choices are qualitatively better than a truncated 
asymptotic expansion. An intuitive explanation is that the rational function can have poles at the 
same points that the function to be approximated does. Because of this, the range of validity of 
the expansion may be extended. 

In a simple version of resummation, we might try to represent only a few terms in JD^ up<J , say 
out to order e 2 , as a rational fraction of operators. The first few terms from (122) are 



jjjpupa _ T livpo 

The [0, 1] order Pade approximation is 



e 2 



eu a d a + (eu a d a ) 2 - -h^d a p + ■■■ 



(125) 



T pupcr 

jQiwpo = n 26 ) 

1 + eu a d a y ' 



while the order [0, 2] approximation can be expressed as 



T 



pupa 



jjjpvpa = ; / 127 x 

[1 + eu<*d a + ±e 2 h^d KX } V ; 



These representations of the viscous stress tensor should already provide a considerable im- 
provement on the standard diffusive form, but we would like to suggest what we think is yet a 
better version. In the earlier work on the thermal smoothing problem, it was possible to obtain an 
operator that had very accurate limits at large and small e. In the normal case, where the transition 
between the two situations is smooth, this meant that the approximation would be reasonably good 
for all values of e. We would like to achieve a similar accuracy for viscous smoothing and, since 
(127) is always good for vanishing e, we need to focus on large e in our determination. For this, we 
propose the following generalization of the prescription for fixing the numerical coefficients in the 
rational approximation. 

We introduce a new rational approximation 



(128) 
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where, as before, P m (e) and Q n (e) are polynomials in e whose orders are to be specified at our 
convenience. However, this time, the coefficients are to be chosen according to a more general 
matching prescription. We suppose that the asymptotic series for 9?^ n ] is 



i=0 



and we determine its coefficients by the conditions that R\ N ^ = 7« for i = 0, 1, , ...,m — 1 and for 
i = N + m,...,N + n + m — l where the orders m and n are to be specified as before. Thus we leave 
a gap between the highest and lowest order terms that we choose to match with the asymptotic 
series of the original function. 

By varying N we control how much influence the higher terms in the series can have on our 
rational approximation. On the other hand, nothing has been lost since it is clear that when 
N = this procedure reduces to the standard one: , = K[ m>n ] . However, we shall be especially 
interested in 

d co \[g( e )}= lim d N) Me)] . (130) 
For example, if we represent g(e) at the level of approximation where we have 

<ilb(<>] = 7-^ (131) 



we obtain 



and 



We have then, in particular, that 



with 



Since lim (7o) 1/W = 1 for any finite 70, we have 



/ \ 1 / N 

7AT+1 nr>0 , 

q N = . (133) 

7JV 

1 yoo fc 

l/N 

lim ( -iil ) . (135) 

iV->oo V 70 / 



N^oo 



qQO = lim ( 7JV )V" = lim 2^±i . (136) 

N— >oo Af— >oo 

For the present study, we shall use only the representation for the viscous smoothing operator 
with N = 00 and [m, n] = [0, 1]. It is therefore convenient to use a simple notation for this case, 
namely 

=<?][•••]■ (137) 
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Returning now to the study of the viscous smoothing operator in (122), we apply the resum- 
mation procedure to find that 



K 



1 + eu a d a 



(138) 



.71=0 

The disappearance of j from this expression is a result of the large N limit. 



Then we note that 

oo 

g^-i^-ITi^- (139> 

In terms of these newly defined smoothing operators, we may write 

jjjfiupa = t W VxVt ( 140 ) 

where ^ 

Vx = i + le*h"Pd af3 ' Vt = l + eu-8 a • 

Then the viscous stress tensor can be written as 

3"" = ^t^I^Si^) (142) 

or, more compactly, within our present level of approximation, as 

E^ = V x V T [fi r^u p , a ] (143) 

where hq = 16 ^ S is the first order shear viscosity of the Thomas theory 

If we denote the Thomas stress tensor of (33) by 

= MT^ pa u p , a , (144) 

then, from (143), (138) and (139), we arrive at this linear partial differential equation for E flu : 

(1 + I e 2/ l a/3^)(l + eu " da )Z^ = Eg . (145) 

This is an evolution equation or, what meteorologists would call, a prognostic equation. However, 
the operator u^d^ is effectively a time derivative and, for times longer than the flight times of 
photons, it is negligible in this formula. So when we are dealing with nonrelativistic motions, (145) 
reduces to the diagnostic equation 

^h^dapET = 3£f - E^ . (146) 
An illustration of an application of this formula is given elsewhere (Chen PhD thesis). 
In the nonrelativistic case, (146) could be written as 

-e 2 AE^ u = E£? - E^ (147) 
3 

where A is the three-dimensional Laplacian operator. However, this last step involves the neglect 
of retardation times, which is not always a safe practice (Delache and Froeschle 1972). 
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11. Transport Coefficients 

Our formula for the radiative stress tensor provides a closure relation that may be used in 
the standard evolution equations for the radiative energy density and flux. These equations, in 
conjunction with the fluid dynamical equations, provide a two-fluid description of radiative fluid 
dynamics (Simon f963; Hsieh and Spiegel 1976) that is reasonably complete. However, the simpler 
one-fluid formulation is more tractable and could be used in certain cases, as when radiative forces 
are weak. 

In the single fluid approximation we may think of the radiation as providing a source of 
dissipation for the fluid. In that case, we would treat the fluid as imperfect with appropriate 
dissipative terms and use the transfer theory to provide the corresponding transport coefficients. 
This version of radiative fluid dynamics has been developed by Weinberg (1972). We add to this 
theory here by providing some resummed results for the coefficients of viscosity and conductivity. 

Since the procedures we have laid down provide formal solutions to the transfer equation, we 
may compute the radiative quantities in terms of the material properties of the radiating medium. 
We have described the calculation of the viscous stress tensor since that is the object needed to 
close the hierarchy of moment equations at a low order. We can similarly solve for the energy 
density and radiative flux. 

On resumming the expansions for the energy density and flux of the radiation we obtain 
expressions for them like that for the viscosity tensor. We shall not give the calculations for these 
other moments since they parallel those given above and, in any case, our intention in this section 
is merely to indicate this other possible direction for treating the dynamics of radiating fluids. So 
we simply report that 

A-TT 

E = —V x V T [(S-4euP p )S] (148) 

and 

= -^V X V T [h™ (S p - ASu a u p , ff )\ . (149) 
These expressions may be introduced into the radiative stress tensor, which we write here as 

= Tg v + AT^ (150) 

where the equilibrium part is given by 

Air 

Tq U = AirSu^u" - -^-Shr (151) 
and the nonequilibrium contribution is 

AT"" = {E- AttS^u" + F»u v + F»u» - \{E - 4xS)br + H* 4 " . (152) 
With the formula for E just given, we find that 

E — AirS = ^^V x V T {WS, a + ASv? p ) , (153) 
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while F" is given by (149). 

The transport coefficients can be identified from the expression for entropy generation. This 
is given by Weinberg (1972) as 

S* = AT"" - ^) (154) 

where T is temperature, the entropy current is 

= su» - T- VAT"" (155) 

and s is the total (matter plus radiation) entropy density in a frame locally moving with the matter, 
as advocated by Mihalas (1984). (We have made some sign modifications from Weinberg to adapt 
these expressions to the metric signature used here.) 

The calculations needed in the evaluation of the entropy generation are not difficult, but they 
are lengthy and not central to this work, so we omit them here and merely quote the result, which 
is 



2u 2 2 

+ -fh^h^iu^ + - -V»)(V + U °,P ~ 3 h t» e ) ( 156 ) 



where 6 = and the transport coefficients are given by 



^ ^ ( 167re£\ ,„ 
r] = V X V T —— , (157) 



V X V T 



V 3T J ' 

167reS" fu a T a 1 

— + - 

TO 3 



(158) 



H = V X V T ( ^eS ) . (159) 



and 

/;, = T>JD„ ( 

The positivity of the three generalized transport coefficients is a sufficient condition for the entropy 
generation rate to be positive. This is of interest in that a failure of this rate to be positive would 
signal a breakdown in the procedures. 

The first term in (156) represents entropy generation by thermal processes — heat flux, mo- 
mentum exchange with the radiation field — and the corresponding transport coefficient, r], is a 
generalized thermal conductivity. The second term comes from entropy generation due to the bulk 
volume expansion, and its coefficient, £, is similarly a form of bulk viscosity. Finally, the third term 
results from dissipation by shearing motions, with fi as a generalization of the Thomas formula. 
Each of these transport coefficients take the form of a bare coefficient (the coefficient obtained in 
the diffusive limit) renormalized by the operation of V X V T . However, since these formulae are to 
be used in connection with a stress tensor of the diffusive form, these results are applicable only 
when the flow does not oscillate rapidly in time or space. 
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12. Conclusion 

In the moment method of radiative transfer (Thorne 1981; Struchtrup 1997), one replaces the 
transfer equation by an infinite hierarchy of moment equations. To render this system usable, 
one truncates it by introducing a supplementary closure relation among the retained moments. 
That approximation is typically obtained by truncating the expansion for the intensity at an order 
corresponding to the level of closure desired. When we do this at the lowest reasonable level, only 
the first two moment equations are retained, and these need to be closed by an expression for the 
pressure tensor. The standard procedure at this level, the Eddington approximation, fails badly at 
representing viscous effects in radiative fluid dynamics (Anderson and Spiegel 1972). 

Here, we are proposing to keep the first two moment equations since these contain the main 
physics of the radiative action on the fluid. But rather than following the procedure of truncating 
the moment hierarchy with the corresponding Eddington approximation, we seek a closure approx- 
imation by another route. Our procedure is to develop an asymptotic expansion for the intensity in 
the manner of the Hilbert expansion of kinetic theory (Grad 1963; Caflisch 1983), but we use this 
solution only for the purpose of devising an expression for the viscosity tensor; the detailed solution 
itself does not play a role in the theory. Since the Hilbert expansion is an outer expansion, in the 
sense of asymptotic theory (Bender 1978), we extend its validity by the technique of resummation 
so as to render it useful for situations where the photon mean free path is long 

Because of the structure of the transfer equation, our series for the distribution function turns 
into a development of en^d^, the operator that appears in the transfer equation (20) and that 
operates iteratively on the source function in the expansions. To extract the desired information 
from this development we needed to perform integrals over the momentum space of the photons 
and this required extensive manipulations. Since the expansion parameter e appears in front of 
the highest (indeed the only) derivatives in the problem, we were faced with a problem in singular 
perturbation theory, a subject where successful approximation techniques are known but whose 
rigor is rarely, if ever, established. 

Even with the help of the asymptotic methods, and resummation, the calculations were de- 
manding and we were led to idealize the problem to the case with constant mean free path, though 
this approximation could be avoided by a suitable choice of optical coordinates, as we noted. The 
resummations, in conjunction with approximations, permitted us to obtain a compact expression 
for the radiative viscosity tensor for use at for all photon free paths. Our form for this tensor can 
be written as an integral operator on the diffusion tensor of Thomas (1930). Since this integral 
operator is simple, we could see at once how to convert the result into the linear partial differen- 
tial equation (145) for E^ u . Depending on whether the fluid dynamics is very relativistic or not, 
this equation can be either prognostic or diagnostic but, in any case, it provides a general closure 
relation for radiative fluid dynamics. 

Asymptotic procedures, when suitably applied, can give very good results, and they have 
proved their usefulness in kinetic theory (Caflisch 1983). But of course the results need to be 



-28- 



tested, for example, by comparison with the results of other methods. To do this, we must solve 
specific problems with the methods being examined and to compare results with those from other 
methods. A first test of this approach would be to compare it with the moment method for this 
problem since that theory has been under serious development (Thorne 1981; Struchtrup 1997) and 
should provide a good benchmark for new procedures. We know already that to get anything like 
a reasonable representation of viscous effects in a moment method, we need to truncate at a higher 
order than in the standard Eddington approximation (Anderson and Spiegel 1972). 

Of course, both these methods give macroscopic descriptions of radiative fluid dynamics and 
so cannot satisfy the full boundary and initial conditions of transport theory, as is known in kinetic 
theory (Grad 1963). But we suggest that this loss is more than compensated by a gain in flexibility 
and usefulness of our macroscopic description. This we have seen in applying our method to the 
problem of entropy generation in cosmology (to be reported elsewhere). Those calculations involve 
relatively smooth flows and are not as demanding of a method as more complicated radiating flows 
such as may arise in the dynamics of photon bubbles or radiative interaction with vortices and flux 
tubes in hot media. We believe that the study of such more complicated flows involving strong, 
possibly turbulent, fluctuations in radiating media will provide the real test of our approach. 

Finally, we would suggest that the methods we have used here could be useful in other problems 
where the mean free path of the basic elements is long compared to scales of variation in the 
medium. An example is turbulent convection where the local mixing length theory encounters a 
loss of validity because the mixing length may be comparable to scales of variation in space and 
time in the medium. The results derived here may point the way to a more useful theory of eddy 
transport processes. 



When we introduce the angular expansions into the integrals for the moments of the intensity 
we naturally encounter a variety of integrals over products of the P and nP 1 alone as well as over 
various mixed products of them. In this appendix we list some formulas for such products that we 
have had need of throughout this work. Such formulas may be obtained by writing them in terms 
of linear combinations of u^u v , rf* , and so on, and using special cases to find the coefficients as in 
previous studies (Anderson and Spiegel 1972; Thorne 1981; Struchtrup 1997, 1998). 

First we note that the angular integral of an odd number of factors vanishes. Therefore, for 
the general moment over factors of V 1 we consider (81) only for even £. As can be seen, M^ 1 -^ is 
orthogonal to in all suffices and this simplifies the problem. 

The special cases for £ = 1 and 1 = 2 are most commonly encountered and so it is worth 
writing them out explicitly: 



A. Useful Angular Integrals 




(Al) 
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[ m v in a <m = —{w v hr + w p h vc + hrh?p) . (A2) 

J 15 

The general case is found by induction and we have 

Afr 

M m-iMu = h Qni» . . . h iMu-ii*u) : (A3) 

2^ "I - 1 

where the symmetrization is assumed for the indices inside the parentheses. 
We also make use of moments like 

the simplest of which is 

J n»n u dn = ^-(4u^u u - if) = y (3u^u u - h» v ). (A5) 
The general case is easily obtained using the binomial theorem. On noting that 

= J2(l) n(M1 • • • u ^ l ^ +l ■ ■ ■ l ^ k) ■ ( A6 ) 



i=0 



1 1 fc! 



where binomial coefficient is I ^ 1 = j^izfj] we h ave 

i=o W 

B. The Operators X and £ 

In Section 9 of the text, we have used the following approximation 

oo oo oo oo 

*> = E E ^,™}(-^) m+2 ^{^} = E E c M (-er^c {jM (bi) 

j'=0 m=0 j=0 m=0 

where C{_j,m} represents the rest coefficient in operator ID. To justify this approximation, let us 
define the difference between TL and C as 

L{j,m} = lL{j,m} - £{j,m} ■ (B2) 

Evidently (as we saw in the text), £{o,o} = -^{1,0} = ^{0,1} = 0- This leads to 

OO OO OO CO 

J2J2(- e ) m+2jC {j,m} L {j,m} = ^^(-e) m+23 % ra }i{ 3 , m } 
j=0 m=0 j=2 m=0 

00 

+ ^2 (~ e ) mC {0,m} L {0M 
m=2 

oo 

+ J2(-ef +2 C {lM L {lM . (B3) 

m=0 



After shifting indices as we did previously, we get 

oo oo oo oo 

E ^l2(- e ) m+23C {jM L {j,m} = E E(~ £ ) m+2J ,( ? C {j+l,m} L {j+lrt 
j=0 m=0 j=l m=0 

oo 

+ E (~ £ ) me2 ( C '{0,m+2}^{0,m+2} + C{l,m}^{l,m} ) • ( B4 ) 
m=0 

We can then decompose the summation for £{j, m } as 

oo oo oo oo oo 

E E (~^ m+23C {JM C {jM = EE (- e ) m+23C {j,m} C {JM + E (- e )" 1C '{0,m}' C {0,m})- (B5) 
j'=0 m=0 j=l m=0 m=0 

Thus with the decomposition of (B4)-(B5) we have 

oo oo 

*> = EE%w(- e ) m+2j [% m} + % m }] 

j=Q m=0 

oo oo 

= EE(- £ r +2,c M( £ M+ f2i M) ( B6 ) 

j=l m=0 

oo 

+ E(~ [^{0,m}^{0,m} + e2 (C{0,m+2}^{0,m+2} + C{l,m}^{l,m})] • 
m=0 

We see that the term by term difference between JL and C is of order e 2 . If we omit in each 
term the higher correction associated with L, we are left with 

oo oo oo 

10 = ^(-e) m+2, Cj j)m) £j j)ffl) I ^(-f) m C{o im) i;{o )m ) 

j=l m=0 m=0 

oo oo 

= EE(- f ) m+2Jc w%™)- ( B7 ) 

This is the approximation we used in the text. 
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